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ABSTRACT 

We analyze the mixed frame equations of radiation hydrodynamics under the approximations of flux- 
limited diffusion and a thermal radiation field, and derive the minimal set of evolution equations that 
includes all terms that are of leading order in any regime of non-relativistic radiation hydrodynamics. 
Our equations are accurate to first order in v/c in the static diffusion regime. In contrast, we show 
that previous lower order derivations of these equations omit leading terms in at least some regimes. 
In comparison to comoving frame formulations of radiation hydrodynamics, our equations have the 
advantage that they manifestly conserve total energy, making them very well-suited to numerical 
simulations, particularly with adaptive meshes. For systems in the static diffusion regime, our analysis 
also suggests an algorithm that is both simpler and faster than earlier comoving frame methods. We 
implement this algorithm in the Orion adaptive mesh refinement code, and show that it performs well 
in a range of test problems. 

Subject headings: hydrodynamics — methods: numerical — radiative transfer 



1. INTRODUCTION 

Astrophysical systems described by radiation hy- 
drodynamics span a tremendous range of scales 
and parameter regimes, from th e interiors of stars 
(e.g. iKippenhahn fc WeigertI 1 19941) to accre t ion d isks 
around compact objects (e.g. iTurner et all l2003f l to 
dusty accretion flows arou nd massive protostars (e.g. 
iKrumholz et al.l [20051120071) to g alactic-scale flows onto 
AGN (e.g. [Thompson et al.ll2005[ ). AU of these systems 
have in common that matter and radiation are strongly 
interacting, and that the energy and momentum carried 
by the radiation fleld is significant in comparison to that 
carried by the gas. Thus an accurate treatment of the 
problem must include analysis of both the matter and 
the radiation, and of their interaction. 

Numerical methods exist to simulate such systems in a 
variety of dimensionalities and levels of approximation. 
In three dimensions, treatments of the matter and radi- 
ation fields generally adopt the flux-limite d diffusion ap- 
proximation, flrst introduced bv lAlme fc W ilson (1973), 
for reasons of co mputational cost and simplicity (e.g. 
iHaves et al.ll2006[) . Flux-limited diffusion is optimal for 
treating continuum transfer in a system such as an accre- 
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tion disk, stellar atmosphere, or opaque interstellar gas 
cloud where the majority of the interesting behavior oc- 
curs in optically thick regions that are well described by 
pure radiation diffusion, but there is a surface of optical 
depth unity from which energy is radiated away. Apply- 
ing pure diffusion to these problems would lead to un- 
physically fast radiation from this surface, so fiux-limited 
diffusion provides a compromise that yields a computa- 
tionally simple and accurate description of the interior, 
while al so giving a rea sonably accurate loss rate from the 
surface ()Castoiil200l) . 

However, the level of accuracy provided by this ap- 
proximation has been unclear because the equations of 
radiation hydrodynamics for flux-limited diffusion have 
previously only been analyzed to zeroth order in v/c. 
In contrast, several authors have analyzed the radiation 
hydrodynamic equatio ns in the general case to beyond 
first order in v/c (e.g. iMihalas fc Weibel-Mihalai 119991 : 
ICastoij |2004 and references therein) . In a zeroth or- 
der treatment, one neglects differences between quanti- 
ties in the laboratory frame and the comoving frame. 
The problem with this approach is that in an optically 
thick fluid, the radiation flux only follows Fick's law 
(F cx —VE) in the comoving frame, and in other frames 
there is an added adv ective flux of r adiation enthalpy, as 
first demonstrated bv ICastorl (|1972D . In certain regimes 
(i.e. the dynamic diffusion limit - see b elow) this advec- 
tive flux can dominate the diffusive flux (|Mihalas fc Aueil 
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does give a flux-limiter usable to 
first order in w/c, which is an approach to the problem of 
flux-limiting with relativistic corrections that is an alter- 
native to the one we pursue in this paper. However, this 
approach does not correctly handle the dynamic diffusion 
limit, a case that as we show requires special attention 
becau se order l (? terms can be important. Further- 
more, iPomraningj derives his flux-limiter directly from 
the transfer equation, so the computation provides little 
insight into the relative importance of radiation hydrody- 
namic terms, and the level of accuracy obtained by using 
the uncorrected flux-limiter, the most common procedure 
in astrophysical ap plications. 

iMihalas h KleinI (1982) were the first to derive 
the mixed frame equations of radiation hydrodynam- 
ics dynamics to order vjc in frequency- integrated 
and frequency-dependent forms, and gave n u merica l 
algorithms for solvin g the m. iLowrie et al.l (fl999l) , 
iLowrie fc Morell (|2001h . and 'Huben v fc Burrows! ()2006l ) 
give alternate forms of these equations, as well as numer- 
ical algorithms for solving them. However, these treat- 
ments require that one solve the radiation momentum 
equation (and for the frequency-dependent equations cal- 
culate over many frequencies as well), rather than adopt 
the flux-limited diffusion approximation. While this 
preferable from a standpoint of accuracy, since it allows 
explicit conservation of both momentum and energy and 
captures the angular-dependence of the radiation field in 
a way that diffusion methods cannot, treating the radia- 
tion momentum equation is significantly more computa- 
tionally costly than using flux-limited diffusion, making 
it difhcult to use in three-dimensional calculations. 

In this Paper we analyze the equations of radiation 
hydrodynamics under the approximations that the radi- 
ation field has a thermal spectrum and obeys the flux- 
limited diffusion approximation, and that scattering is 
negligible for the system. Our goal is to derive an ac- 
curate set of mixed frame equations, meaning that radi- 
ation quantities are written in the lab frame, but fluid 
quantities, in particular fluid opacities, are evaluated in 
the frame comoving with the fluid. This formulation is 
optimal for three-dimensional simulations, because writ- 
ing radiation quantities in the lab frame lets us use an 
Eulerian grid on which the radiative transfer problem 
may be solved by any number of standard methods, while 
avoiding the need to model the direction- and velocity- 
dependence of the lab frame opacity and emissivity of a 
moving fluid. 

In § [2] we begin from the general lab frame equations 
of hydrodynamics to first order in vjc, apply the flux- 
limited diffusion approximation in the frame comoving 
with the gas where it is applicable, and transform the ap- 
propriate radiation quantities into the lab frame, thereby 
deriving the corresponding mixed frame equations suit- 
able for implementation in numerical simulations. We re- 
tain enough terms to ensure that we achieve order unity 
accuracy in all regimes, and order vjc accuracy for static 
diffusion problems. In § [3] we assess the significance of 
the higher order terms that appear in our equations, and 
consider where treatments omitting them are acceptable, 
and where they are likely to fail. We show that, in at 
least some regimes, the zeroth order treatments most 
often used are likely to produce results that are incor- 



rect at order unity. We also compare our equations to 
the comoving frame equations commonly used in other 
codes. In § 2] we take advantage of the ordering of terms 
we derive for the static diffusion regime to construct a 
radiation hydrodynamic simulation algorithm for static 
diffusion problems that is simpler and faster than those 
now in use, which we implement in the Orion adaptive 
mesh refinement code. In § [5] we demonstrate it in a 
selection of test problems. Finally, we summarize our 
results in § ini 

2. DERIVATION OF THE EQUATIONS 

In the discussion that follows, we adopt the conven- 
tion of writing quantities measured in the frame comov- 
ing with a fluid with a subscript zero. Quantities in 
the lab frame are written without subscripts. We write 
scalars in italics (e.g. a), vectors in boldface (e.g. a), 
and rank two tensors in calligraphy (e.g. A). We indi- 
cate tensor contractions over a single index by dots (e.g. 
a • b = a*&i), tensor contractions over two indices by 
colons (e.g. A:B — A^^ Bij), and tensor products of vec- 
tors without any operator symbol (e.g. (ab)*-' — a^V). 
Also note that we follow the standard convention in ra- 
diation hydrodynamics rather than the standard in astro- 
physics, in that when we refer to an opacity k we mean 
the total opacity, measured in units of inverse length, 
rather than the specific opacity, measured in units of 
length squared divided by mass. Since we are neglecting 
scattering, we may set the extinction x = k. 

2.1. Regimes of Radiation Hydrodynamics 

Before beginning our analysis, it is helpful to examine 
some characteristic dimensionless numbers for a radia- 
tion hydrodynamic system, since evaluating these quan- 
tities provides a useful guide to how we should analyze 
our equations. Let £ be the characteristic size of the 
system under consideration, u be the characteristic ve- 
locity in this system, and Ap ^ 1/k, be the photo n mean 
free path. Following [Mihalas fc Weibel-Mihalai (fl999). 
wc can define three distinct limiting cases by considering 
the dimensionless ratios r = €/Ap, which characterizes 
the optical depth of the system, and (3 = u/c, which 
characterizes how relativistic it is. Since we focus on 
non- relativistic systems, we assume /3 <C 1. We term the 
case T <C 1, in which the radiation and gas are weakly 
coupled, the streaming limit. If r ^ 1 then radiation and 
gas are strongly coupled, and the system is in the diffu- 
sion limit. We can further subdivide the diffusion limit 
into the cases /3 ^ and (3 <C r^^. The former is 
the dynamic diffusion limit, while the latter is the static 
diffusion limit. In summary, the limiting cases are 

T ^ 1, (streaming limit) (1) 

T > 1, < 1, (static diffusion hmit) (2) 
T 3> 1, /?T 3> 1, (dynamic diffusion limit). (3) 

Physically, the distinction between static and dynamic 
diffusion is that in dynamic diffusion radiation is prin- 
cipally transported by advection by gas, so that terms 
describing the work done by radiation on gas and the 
advection of radiation enthalpy dominate over terms de- 
scribing either diffusion or emission and absorption. In 
the static diffusion limit the opposite holds. A paradig- 
matic example of a dynamic diffusion system is a stellar 
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interior. The optical depth from the core to the surface 
of the Sun is r ~ 10^^, and typical convective and rota- 
tional velocities are ^ 10~^^c = 0.3 cm s~^, so the Sun 
is strongly in the dynamic diffusion regime. In contrast, 
an example of a system in the static diffusion limit is a 
relatively cool, dusty, outer accretion dis k around a form- 
ing m assive protostar, as studied e.g. bv lKrumholz et all 
()2007( ). The specific opacity of gas with the standard in- 
terstellar dust abundance to infrared photons is n/ p ^ 1 
cm^ g~^, and at distances of more than a few AU from 
the central star the density is generally p ^ 10""'^^ g cm~^. 
For a disk of scale height h ^ 1{) AU, the optical depth 
to escape is 

-1 



k/p 



2 —1 

cm g 



]^Q-i2 g cm ^ 



h 



10 AU 



The velocity is roughly the Keplerian speed, so 



/3« 1.4 X 10" 



1/2 



-1/2 



(4) 



(5) 



IQMqJ VlOAU/ 

where is the mass of the star and r is the distance 
from it. Thus, this system is in a static diffusion regime 
by roughly two orders of magnitude. 

In the analysis that follows, our goal will be to obtain 
expressions that are accurate for the leading terms in all 
regimes. This is somewhat tricky, particularly for diffu- 
sion problems, because we are attempting to expand our 
equations simultaneously in the two small parameters [3 
and 1 /r. The most common approach in radiation hydro- 
dynamics is to expand expressions in powers of (3 alone, 
and only analyze the equations in terms of r after drop- 
ping terms of high order in (3. However, this approach 
can produce significant errors, because terms in the radi- 
ation hydrodynamic equations proportional to the opac- 
ity are multiplied by a quantity of order r. Thus, in our 
derivation we will repeatedly encounter expressions pro- 
portional to /3^T, and in a problem that is either in the 
dynamic diffusion limit or close to it {(3t ^1), it is in- 
consistent to drop these terms while retaining ones that 
are of order (3. We therefore retain all terms up to order 
0^ in our derivation unless we explicitly check that they 
are not multiplied by terms of order t, and can therefore 
be dropped safely. 



2.2. The Equations of radiation hydrodynamics 

We now start our derivation, beginning from 
the lab frame equations of radiation hydrodynamics 
Mihalas fc Klein" 1 9821 : iMihalas fc Weibel-Mihalaslfl999l : 
Mihala s & Auer. .200lF 

dp 



dt 



4- V • {p^r) = 



d_ 
dt 



(PV) + V • (PVV) : 



-VP 



d_ 

dt 



{pe) + W ■[{pe + P)v]^cG" 



dE 
'dt 

1 ar 

'^~dt 



V • F: 



G 



(6) 
(7) 
(8) 
(9) 
(10) 



where p, v, e, and P are the density, velocity, specific 
energy (thermal plus kinetic), and thermal pressure of 
the gas, E, F, and V are the radiation energy density, 
flux, and pressure tensor. 



cE^ dv I dni{n,i^ 



dv i dfi n/(n, v) 



(11) 
(12) 



cP = j dv j dnnnl{n,u), (13) 
(G'^, G) is the radiation four- force density 



CG°: 



cG 



dv / dn [K(n, v)I{n, v) - r]in, v)], (14) 



dv / dfl[K{n,v)I{n,v) — ri{n,v)]n, (15) 



and /(n, v) is the intensity of the radiation field at fre- 
quency V traveling in direction n. Here ^(n, v) and 
77(11, v) are the direction- and frequency-dependent ra- 
diation absorption and emission coefficients in the lab 
frame. Intuitively, we can understand cG*^ as the rate 
of energy absorption from the radiation field minus the 
rate of energy emission for the fluid, and G as the rate 
of momentum absorption from the radiation field minus 
the rate of momentum emission. Equations ([SI) - ([8]) are 
accurate to first order in v/c, while equations ^ - (jlOp 
are exact. Note that no terms involving opacity or opti- 
cal depth appear explicitly in any of these equations, so 
the fact that they are accurate to first order in f3 means 
that they include all the leading order terms. 

In order to derive the mixed- frame equations, we must 
now evaluate the radiation four- force (G°,G) in terms 
of lab frame radiation quantities and comoving frame 
emission and absorption coefficients. IMihalas fc AueH 
([2(301) show that, if the flux spectrum of the radiation 
is direction-independent, the radiation four- force on a 
thermally-emitting material to all orders in v/c is given 
in terms of moments of the radiation field by 

G^^lh^KQE + (1 - j'^)kqf]E - jKQpanTQ 

- 7(v • F/c^)[kof - 27^(kof - kqe)] 
-73(kof -koe)(vv):7'/c2, (16) 

G = 7Kof(F/c) - 1 kqp a rT^ {-v/c) 

- [7^(«0F - koe)(v/c)£; + 7Kof(v/c) • V] 

+ j''{kof - «oe)[2v • F/c^ - (vv):7'/c3]v,(17) 



where 7 = 1/ -\/l — w^/c^ is the Lorentz factor and Tq 
is the gas temperature. The three opacities that appear 
are the Planck-, energy-, and flux-mean opacities, which 
are defined by 

dvo Ko{vo)B{vo,To) 



Kqp = ■ 



KOF : 



Bin) 

dvQ Ko{vo)Eo{vo) 
E'o 

dvo Ko{vo)Fo{vo) 



(18) 



(19) 



(20) 



where Eq{vq) and Fo(t'o) are the comoving frame ra- 
diation energy and flux per unit frequency, Eq and Fq 
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are the corresponding frequency-integrated energy and 
flux, and B{iy,T) = {2hiy^ /c^)/{e'"'^''^^ -I) and B(T) = 
cauT'^ / (An) are the frequency-dependent and frequency- 
integrated Planck functions. 

Note that we have implicitly assumed that the opacity 
and emissivity are directionally-independent in the fluid 
rest frame, which is the case for any conventional ma- 
terial. We have also assumed that the flux spectrum is 
independent of direction, allowing us to replace the flux- 
mean opacity vector with a scalar. This may not be the 
case for an optically thin system, or one in which line 
transport is important, but since we are limiting our ap- 
plication to systems to which we can reasonably apply 
the diffusion approximation, this is not a major limita- 
tion. 

To simplify {G''\ G), first we assume that the radiation 
has a blackbody spectrum, so that £^0(1^0) oc B(vq,Tq). 
In this case, clearly 

KOE = Kop. (21) 

Second, we adopt the flux-limited diffusion approxima- 
tion (see below), so in optically thick parts of the flow 
Fo(i^o) — V£'o(i^o)/Ko(i^o) (Fick's Law). This implies 
that Fo(i'o) « -[dB{iyo,To)/dTo]{VTo)/KoM, and sub- 
stituting this into ([20|) shows that the flux-mean opacity 
kqf is equal to the Rosseland-mean opacity, defined by 



" V V 




-E+--V 


4( 


. c c . 





-1 



Jq diyQKoiiyo) — ^g^^ 



3('^0:To) 
dTo 



(22) 



In optically thin parts of the flow, |Fo(i'o)| cEq^vq), 
so in principle we should have kqf = kqe- However, in- 
terpolating between these cases is complex, and the flux- 
limited diffusion approximation is of limited accuracy for 
optically thin flows in complex geometries. Moreover, 
our approximation that the radiation spectrum is that 
of a blackbody at the local radiation temperature is it- 
self problematic in the optically thin limit, so setting 
kqf = kqp would not necessarily be more accurate than 
using kqr- We therefore choose to optimize our accuracy 
in the optically thick part of the flow and set 

kqf = KOR- (23) 

With these two approximations, the only two opacities 
remaining in our equations are kqr and kqp, both of 
which are independent of the spectrum of the radia- 
tion field and the direction of radiation propagation, and 
which may therefore be tabulated as a function of tem- 
perature for a given material once and for all. 

Next, we expand (G'',G) in powers of v/c, retaining 
terms to order jf?. In performing this expansion, we 
note that |F| < cE, and Tr(-p) = E. The resulting 
expression for the radiation four-force is 

= Kop [e - ) + (kor - 2kop) 



1 iv\^ 

2 \c) 



2(kop - kor)-B + Kop LE - 



47rB 



vv / 

-I- (kqp - kor)-^:?^ + 0\-^ 



(24) 



G = KoR— -I- Kop (^-J \ E — 1 



2(kor - Kop 



(v • F)v 



(25) 



It is helpful at this point, before we making any fur- 
ther approximations, to examine the scalings of these 
terms with the help of our dimensionless parameters /? 
and T. In the streaming limit, radiation travels freely 
at c and emission and absorption of radiation by matter 
need not balanc e, so |F| ^ c-E and 47ri3/c — E ~ E . For 



static diffusion, iMihalas fc Weibel-MihalasI (|1999D show 
that |F| - cEjr and 47rB/c - E E/t'^. For dy- 
namic diffusion, radiation travels primarily by advection, 
so |F| ^ vE. We show in Appendix [X] that for dy- 
namic diffusion AttB jc — E^ (?E. Note that the scaling 
AttB/ c-E ~ {f3/T)E given in IMihalas fc Weibel-Mihalad 
(|1999( ) appears to be incorrect, as we show in the Ap- 
pendix. Using these values, we obtain the scalings shown 
in Tabled] for the terms in ^ and ([25]) . 

The Table shows that, despite the fact that we have 
kept all terms that are formally order or more, in 
fact we only have leading-order accuracy in the dynamic 
diffusion limit, because in this limit the order unity and 
order (3 terms in G° vanish to order (3'^. To obtain the 
next-order terms, we would have had to write G° to order 
f3^. A corollary of this is that treatments of the dynamic 
diffusion limit that do not retain order terms are likely 
to produce equations that are incorrect at order unity, 
since they will have dropped terms that are of the same 
order as the ones that have been retained. 

At this point we could begin dropping terms that are 
insignificant at the order to which we are working, but it 
is cumbersome to construct a table analogous to Table [T] 
at every step of our derivation. It is more convenient to 
continue our analysis retaining all the terms in (|24p and 
(|25|) . and to drop terms only periodically. 

Before moving on, there is a subtlety in ([24]) and (|25p 
that is worth commenting on. Consider a gray fluid, one 
in which kor — kqp = kq. In cG*^, the term that de- 
scribes the work done by radiation, — kov • F/c, has the 
opposite sign from what one might naively expect. Using 
cG° in the gas energy equation ^ in this case implies 
that the gas energy increases when v and F are anti- 
aligned, i.e. when gas moves into an oncoming photon 
flux. We can understand the origin of this somewhat 
counter-intuitive behavior by considering the example of 
a fluid in thermal equilibrium with a radiation field in its 
rest frame (i.e. AttB = cEq). In the comoving frame, the 
radiation four-force behaves as one intuitively expects: at 
leading order the rate at which the radiation field trans- 
fers momentum density to the gas is Go = kqFq/c, and 
the rate at which the gas energy density change s as a 
result is cG[j = koV • Fq/c (|Mihalas fc Aueij|200ll their 
equations 53a and 53b). Thus, gas loses energy when 
it moves opposite the direction of the flux, and hence 
opposite the force. 

However, now consider the fluid as seen by an observer 
in a frame boosted by velocity — v relative to the fluid. 
The observer sees the radiation energy density as E, 
which differs from Eq by 2v • F/c^ (see equation ISTjl. 
and this difference is the reason that the work term in 
G° is — Kov • F/c^. Physically, this happens because an 
observer who sees the fluid moving at velocity v also 
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TABLE 1 

SCALINGS OF TERMS IN THE RADIATION FOUR-FORCE DENSITY 



G" or G 


Term 


Streaming 


Statie Diffusion 


Dy 


namic Diffusion 


G" 


K(,p{E — AivB/c) 


T 


1/x 






GO 


(kor-2kop)(v-F/c2) 


I3t 








G" 


(i;/c)2(kop - kor)-E 


/3V 


ffir 




M ' 


GO 


(l/2)Kc)2«:op(47rS/c-£;) 


/3V 


ffllr 




/3V 


GO 


(kor - kqp){vv / c^y.V 


/3V 


ffir 




/3^x 


G 


kqrF/c 


T 


1 




/3x 


G 


K0p(v/c)(47rB/c-£;) 


f3r 


/3/r 




/33r 


G 


«:or[(v/c)£; + (v/c) -P] 


I3t 


/3r 




/3x 


G 


(1/2)Kc)2«;orF/c 


p-'r 








G 


2{kor - Kop){v ■ F)v/c3 


/3V 


0^ 




/33r 


Note. — Col. (1): Whether the term appears in G" or G. CoL (3)-{5) 
normahzed to E/£. The scaUngs that are of leading order in each regime are 


: All scalings are 
boldfaeed. 



sees the radiation and gas as being out of thermal equi- 
hbrium (47ri? ^ cE), since E and Eq are different. This 
disequihbrium leads the radiation and gas to exchange 
energy at a rate that is opposite in direction and twice 
as large as the radiation work, kqv • F/c. This is why 
the "work" term has the opposite sign than the one we 
might expect. Thus, for the rest of this paper, while for 
convenience we continue to refer to (kqr — 2kop)v • F/c 
and the terms to which it gives rise as "work" terms, it 
is important to keep in mind that in reality this term 
contains contributions from two different effects of com- 
parable magnitude, the "Newtonian work" kqrv • F/c 
and the post-Newtonian term — 2kopv • F/c describing 
the imbalance between emission and absorption that an 
observer sees solely because the fluid is moving. 

With this point understood, we now adopt the flux - 
limited diffusion approximation (|Alme fc WilsonI [19731 ). 
under which we drop the radiation momentum equation 
pU)) and set the radiation flux in the comoving frame to 



cA 



(26) 



where A is a dimensionless number called the flux-limiter. 
Many functional forms for A are possible. For the 
code implementation we descr ibe later, we adopt the 
iLevermore fc Pomranin^ ()1981l ) flux-limiter, given by 



A: 



R = 



— ( coth R - — 
R\ R 

NE, 



01 



'^OR-E'O 



(27) 
(28) 



However, our derivation is independent of this choice. 
Regardless of their exact functional form, all flux lim- 
iters have the property that in an optically thick medium 
A — > 1/3, thereby giving Fq — > — [c/(3KoR)]Vi?o, the cor- 
rect value for diffusion. In an optically thin medium, 
A (KoRi?o/|Vi?o|)no, where no is the unit vector an- 
tiparallel to V£'q, so the flux approaches Fq — > cEono, 
and the propagation speed of radiation is correctly lim- 
ited to c. 

For the lLevermore fc Pomraningl flux-limiter we adopt, 
the corresponding approximate va, lue for the radiation 
pressure tensor is (jLevermord 119841 ) 



Vo 



[(1 - i?2)2: + (3i?2 - l)nono] , (29) 



where X is the identity tensor of rank 2 and 
i?2 = A 



X^R^. 



(30) 



Physically, this approximation interpolates between the 
behavior in very optically thick regions, where i?2 — > 
l/3+0(l/r^), the radiation pressure is isotropic, and off- 
diagonal components vanish, and optically thin regions, 
where i?2 ^ 1 and the radiation pressure tensor is zero 
orthogonal to ng and Eq par allel to it. 

Note tha t for pur e diffusion IMihalas fc Weibel-Mihalai 
(1999) and I Cast oil (12004) show that the pressure tensor 
reduces to {Eo/3)X plus off-diagonal elements of order 
I3/t or Our approximation does not quite reproduce 
this, since in the diffusion limit it gives Po = {Eq/Z)T 
plus off-diagonal elements of order t~^. We might there- 
fore worry that, in the static diffusion regime where 
(3 ^ T~^, we will have an incorrect term. However, 
examination of our final equations below shows that all 
terms arising from off-diagonal elements of Vq are smaller 
than order (3 in the static diffusion limit, so adopting the 
[Levermorc (1984) approximation for the pressure tensor 
does not introduce any incorrect terms at order (3 in the 
final equations. 

To use the approximations ([26]) and (f29l) to evaluate 
the radiation four-force, we must Lorentz transform them 
to express the radiation quantities in the lab frame. The 
Lorentz transforms for the energy, flux, and pressure 
to se cond order \n v/c are (jMihalas fc Weibel-Mihalai 
fl999h 



E = E, 



1 







F = Fo 



vEo - 
vFo 



FqV , 



i;2So + (vv):Po] (31) 
[«2Fo + 3v(vFo)](32) 

v(vPo)]. (33) 



1 
2? 



1 



[vv£'o 



Note that in the expression for V we have simplified the 
final term using the fact that Vo is a symmetric tensor. 

Using the same scaling arguments we used to construct 
Table [Tl we see that V and Vq differ at order f3 in the 
streaming limit, at order (3/t for static diffusion, and at 
order for dynamic diffusion. Since this is below our 
accuracy goal, we need not distinguish V and Vq. The 
same is true of E and Eq. However, F is different. In 
the comoving frame in an optically thick system, one is 
in the static diffusion regime, so Fq ^ cEq/t. Since vEq 
and V • Vq are of order PcEq, and in dynamic diffusion 
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/3 ^ 1/r, this means that vEq and v • Vq are the dom- 
inant components of F in dynamic diffusion, and must 
therefore be retained. Thus, 



cA 



-VE + vE + v - V, 



(34) 



which is simply the rest frame flux plus terms describing 
the advection of radiation enthalpy. 

Substituting with P_= Vq and ^ into the four- 

and continuing to retain 



force density ([24]) and ([251 
terms to order w^/c^, gives 

' AttB 

: KOP I E 



9 

1 /i; 
-XVE- 



c 

3-i?2 



- 2 



A 

3i?2-l 



- 1 



VE 



(v 



KOP- 



E 



E 



AttB 



(35) 



c 



1 /i; 

2 Vc 



AV£; 



2A 



V'^OR 



(v • VE)v 



(36) 



Here n is the unit vector antiparallel to WE. We again 
remind the reader that, although these equations contain 
terms of order /3^, they are not truly accurate to order 
because we did not retain all the when applying 
the Lorentz transform to the flux and pressure. How- 
ever, these equations include all the terms that appear 
at the order of accuracy to which we are working, and 
by retaining terms of order we guarantee that these 
terms will be preserved. 

Inserting (G°, G) and the lab frame flux (f34|) into the 
gas momentum and energy equations ([7]) and ([S]), and 
the radiation energy equation ([9]), and again retaining 
terms to order jc? gives 

d 

-V • (pvv) - VP - WE 



dt 



(PV): 



2A 



- 2 ( c 
(v • \/E)v 



WE 



V «0R / 

[{pe + P)w] - kop{^t:B - cE) 



(37) 



+ A 2 



-l]vVE 



KOP 



E 



3-R2 



3i?2 - 1 



(v-n)2 



^ ) Kop (47rB 



dt 



£; = v 



cE) 

[—\/e] +Kop(47rB 



(38) 



-A 2 



V«OR 
Kqp 



cE) 



KOR 



Kop 



E 



-l\^^■VE 



3i?2 



3-i?2 



-(v.n)2 



3-i?2 



-v£; 



3i?2 



{nn)E 



+ 



i (^)'kop(4^B-c^). (39) 

At this point we construct Table [2] showing the scalings 
of the radiation terms to see which must be retained 
and which are superfluous. In constructing the table, 
we take spatial derivatives to be of characteristic scaling 
i.e. we assume that radiation quantities vary on a 
size scale of the system, rather than over a size scale of 
the photon mean free path. In the streaming limit, A ~ r 
and _R2 ~ 1 + 0(r). In the diffusion limit A ~ 1/3 and 
i?2'- 1/3 + 0(r-2). 

Using Table [5] to drop all terms that are not significant 
at leading order in any regime, we arrive at our final 
equations: 



d_ 
dt 
d_ 
dt 



(pv) = -V • (pvv) - VP - WE 
ipe) = -V • [{pe + F)v] - Kop(47rB 



(40) 



cE) 



+ A 2 



, KOP 
KOR 
i?2 



\IE 



V 

Kqp — E 
c 



(41) 



= V 



VKOR / 



Kop(47rS - cE) 



A ( 2^ - 1 I V • V£; 

KOR 



3 ~ i?2 v\ 

-^-OP-E 
V. {^^vE 



(42) 



These represent the equations of momentum conserva- 
tion for the gas, energy conservation for the gas, and en- 
ergy conservation for the radiation field, which, together 
with the equation of mass conservation ([6]) , fully describe 
the system under the approximations we have adopted. 
They are accurate and consistent to leading order in the 
streaming and dynamic diffusion limits. They are accu- 
rate to first order in /3 in the static diffusion limit, since 
we have had to retain all order (3 terms in this limit 
because they are of leading order in dynamic diffusion 
problems. Also note that if in a given problem one never 
encounters the dynamic diffusion regime, it is possible to 
drop more terms, as we discuss in § |4l 

The equations are easy to understand intuitively. The 
term — WE in the momentum equation (|40p simply rep- 
resents the radiation force kqrF/c, neglecting distinc- 
tions between the comoving and laboratory frames which 
are smaller than leading order in this equation. Similarly, 
the terms ±kop{4:'kB — cE) and ±A(2kop/kor — l)v- Vi? 
in the two energy equations ((4T|) and ((42|) represent radi- 
ation absorbed minus radiation emitted by the gas, and 
the work done by the radiation field as it diffuses through 
the gas. The factor (2kop/kor — 1) arises because the 
term contains contributions both from the Newtonian 
work and from a relativistically-induced mismatch be- 
tween emission and absorption. The term proportional 
to KopE/c represents another relativistic correction to 
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TABLE 2 

SCALINGS OF TERMS IN THE CONSERVATION EQUATIONS 



Equation 


Term 


Streaming 


Static Diffusion 


Dynamic Diffusion 


M 


AVE 


T 


1 


1 


M 


K0p(v/c2)(47rB - cE) 




/3/r 




M 


{l/2){v/c)^X\/E 


fS^T 






M 


2A(«;op/kor-1)(v- V£;)v/c2 


/3V 




/32 


G and R 


Kop{AnB - cE) 


T 


1/t 


/3V 


G and R 


A(2kop/kor - l)v ■ 


fir 


/3 


G and R 


KQp(v^/c)[['i- R2)/2\E 


/3V 




/3^x 


G and R 


Kop[(v-n)2/c][{3i?2-l)/2]i? 


/3V 






G and R 


(l/2){v/cfKop(cE-A-KB) 


/3V 




/3V 


R 


V ■ [{cA/kor)V£;] 


1 


1/x 


1/r 


R 


V-{[(3-iJ2)/2]v£;} 




/3 


/3 


R 


V-{[(3iJ2-l)/2]v(nn)£;} 


13 


/3/r2 





Note . — Col. (1): Which eq uatio n the term appears in. M = momentum (I37I I. G = gas 
energy II38I I. R = radiation energy II39I I. Col. (3)-(5): All scalings are normalized to E/l for the 
momentum equation, and cE/l for the energy equations. Scalings that are of leading order in each 
regime for each equation are boldfaced. 



the work, this one arising from boosting of the flux be- 
tween the lab and comoving frames. In the radiation 
energy equation (|42|). the first term on the left hand side 
is the divergence of the radiation flux, i.e. the rate at 
which radiation diffuses, and the last term on the right 
hand side represents advection of the radiation enthalpy 
E + V hy the gas. 

It is also worth noting that equations (l38|) and (l39|) 
are manifestly energy-conserving, since every term in 
one equation either has an obvious counterpart in the 
other with opposite sign, or is clearly an advection. In 
contrast, the momentum equation (j40p is not manifestly 
momentum-conserving, since there is a force term —AVE' 
with no equal and opposite counterpart. This non- 
conservation of momentum is an inevitable side-effect of 
using the flux-limited diffusion approximation, since this 
approximation amounts to allowing the radiation field to 
transfer momentum to the gas without explicitly tracking 
the momentum of the radiation field and the correspond- 
ing transfer from gas to radiation. 



3. THE IMPORTANCE OF HIGHER ORDER TERMS 

Our dynamical equations result from retaining at least 
some terms that are formally of order 0^ . Even though 
our analysis shows that these terms can be the leading 
ones present, due to cancellations of lower order terms, 
one might legitimately ask whether they are ever phys- 
ically significant. In § 13.11 we address this question by 
comparing our equations to those that result from lower 
order treatments. In § 13.21 we also compare our equations 
with those generally used in comoving frame formulations 
of radiation hydrodynamics. 

To make our work in this section more transparent, 
and since we are more interested in physical intuition 
than rigorous derivation here, we specialize to the diffu- 
sion regime in gray materials. Thus, we set X — R2 = 1/3 
and Kop = kqr = kq. A more general analysis produces 
the same conclusions, but is more mathematically cum- 
bersome. We also focus on the radiation energy equation, 
since all the terms that appear in the gas energy equa- 
tion also appear in it, and because there are no higher 
order terms present in the momentum equation. Under 
these assumptions, our radiation energy equation (j42p 



becomes 
d 



dt ^ (.3^0^ ) 



Ko(47rB - cE) 



4 1 4 

- y.{^E)--^r■VE+-K^-E. (43) 
3 3 3 c 

3.1. Comparison to Lower Order Equations 

A common approach in radiation-hydrodynamic prob- 
lems is to expand the equations in /3, rather than in 
both (3 and r as we have done, and drop at least 
some terms that are of order in every regime (e.g 
iMihalas fc Weibel-Mihalad Il999( ). To determine how 
equations derived in this manner compare to our higher 
order treatment, we compare our simplified energy equa- 
tion (|43p to the corresponding equation one would obtain 
by following this procedure with (|42|) . This resulting en- 
ergy equation is 

l-^E^y\-^VE)+.,iA.B-cE) 

- • {^E) - iv • VE. (44) 

It is important to caution at this point that, as we show 
below, equation |^^[ ) is not accurate to leading order in at 
least some cases, and should not be used for computations 
unless one carefully checks that the missing terms never 
become important in the regime covered by the computa- 
tion. 

Compared to the energy equation that we obtain 
by retaining all leading order terms in /3 and r, (|44p 
is missing the term {A/3)kqv'^E/c. If we think of the 
flux as having two parts, a "diffusion" part proportional 
to Vi? that comes from radiation diffusion in the co- 
moving frame, and a "relativistic" part proportional to 
vE + V ■ V that comes from the Lorentz transformation 
between lab and comoving frames, then it is natural to 
describe the v • WE term as the "diffusion work" aris- 
ing from the combination of the diffusion flux and the 
post-Newtonian emission-absorption mismatch (as dis- 
cussed in § 12. 2p . and the kqv'^E/c as the "relativistic 
work" arising from the relativistic flux. The presence or 
absence of this relativistic work term is the difference be- 
tween our leading order-accurate equation and the equa- 
tion one would derive by dropping terms. Analyzing 
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when, if ever, this term is physically important lets us 
identify in which situations a lower order treatment may 
be inadequate. 

If we use Table [2] to compare the relativistic work 
term to the emission/absorption term, we find that 
{kqv'^E/c)/[kq{A'kB — cE)] is of order /3^r^ for static dif- 
fusion, and of order unity for dynamic diffusion. Thus, 
the term is never important in a static diffusion prob- 
lem, but is always important for a non-uniform, non- 
equilibrium dynamic diffusion problem system. We add 
the caveats about non-uniformity and time-dependence 
because in a system where there is no radiation-gas en- 
ergy exchange, the relativistic work term will be small 
due to a cancellation. The example in Appendix [X] 
shows that in an equilibrium, uniform medium, the terms 
kq{AttB — cE) and (4/3)Ko'f^£'/c cancel exactly at orders 
up to P"^. We expect any system where variations occur 
on a scale for which /3t ^ 1 to resemble such a uniform, 
equilibrium medium, and thus we do not expect the term 
{4: / 3) kqv'^ E / c to be important in such a system. 

That said, there is still clearly a problem with omit- 
ting the relativistic work term in a system where /3t ~ 1. 
In this case. Table [2] implies that every term on the 
right hand side is roughly equally important regard- 
less of whether we use the static of dynamic diffusion 
scalings. To illustrate this point, consider a radiation- 
dominated shock. The width of such a shock is set by 
the balance between radiation diffusing upstream from 
the hot post-shock region into the cold pre-shock region, 
and advection of the radiation back downstream by the 
pre-shock gas. T his condition requires that Pt ^ 1 
across the shock ()Mihalas fc Weibel-MihalasI Il999l ) . so 
its width w ^ Xp/P- Since E changes by of order 
unity across this distance, its spatial derivative is of or- 
der VE ~ E/w ~ {P/\p)E. Applying this to (gS]), we 
find that each term on the right hand side is of order 
P'^{c/Xp)E. Since the terms like -(4/3)V • (vE) describ- 
ing advection and V • [c/(3Ko)Vi?] describing diffusion 
are obviously important in the structure of the shock, 
causing order unity changes in E, and the relativistic 
work term is comparable, it follows that the relativistic 
work term is equally important. One can obtain the cor- 
rect structure within a radiation-dominated shock only 
by retaining the relativistic work term. 

An interesting point to note here is that omitting the 
relativistic work term will not produce errors upstream or 
downstream of a shock, because /3r ^ 1 in these regions. 
Furthermore, the jump conditions across a shock should 
be correct. The omitted term affects radiation-gas energy 
exchange, not total energy conservation, and all that is 
required to get the correct jump conditions are conser- 
vation of mass and energy, plus correct computation of 
the upstream and downstream radiation pressures. The 
lower order treatment will therefore only make errors 
within the shock. Whether this is physically important, 
or it is sufficient to get the jump conditions correct, de- 
pends on whether one is concerned with structures on 
scales for which /3t ~ 1. An astrophysical example of a 
system where one does care about structures on this scale 
is a radiation-dominate d accretion disk su bject to pho- 
ton bubble instability (iTurner et al.l l2003l ) . Such disks 
are in the dynamic diffusion regime over the entire disk, 
but photon bubbles form on small scales within them, 
and individual bubbles may have Pt ^ 1 across them. 



3.2. Comparison to Comoving Frame Formulations 

Many popular numerical treatments of radia- 
tion hydrodynami cs (e. g. iTurner fc Stond 120011 : 

IWhitehouse fc Batd 120041 : iHaves et al.l l2006l ) use a 
comoving formulation of the equations rather than 
our mixed frame formulation. It is therefore useful to 
compare our equations to the standard comoving frame 
equations. In the comoving formulation, the evolution 
equation for the radiation field is usually the first law 
of thermodynamics fo r the comoving radiation field 
(iMihalas & KleinlUflS^ . 

p-^ (y ) + = '*o(4^-S - cSo) - V • Fo. (45) 

This equation is accurate to first order in P in the sense 
that it contains all the correct leading order terms and 
all terms that are smaller than them by order P or less. 

To compare this to our mixed frame radiation energy 
equation ()42|) . we replace the comoving frame energy Eq 
in (|45p with the lab frame energy E using the Lorentz 
transformation (|3ip and retain all terms that are of lead- 
ing order in any regime. In practice, this means that we 
set Eq = E inside the time derivative, since the differ- 
ence between E and Eq is at most /? /r or for static or 
dynamic diffusion. However, when replacing Eq with E 
in the heating/cooling term AttB — cEq, we must retain 
all the terms in (pij) because the leading term AttB — cE 
is itself only of order t^^ or P'^ relative to E, so the dif- 
ference between E and Eq can be of leading order. -'^ This 
gives a transformed equation 

^Wt (7) ^ ^^^^ " «o(4^S - cS) - V • Fo 

+ 2kq^^ + ^ [v^E + (vv): Po] • (46) 
c c 

If we now adopt the diffusion approximation Fq = 
— c/ {3kq)WEq and Vq = {1/3)EqT, use the Lorentz trans- 
formation to replace Eq with E throughout, and again 
only retain terms that are of leading in order in some 
regime, then it is easy to verify that (H^ reduces to 
Thus, our evolution equation is equivalent to the comov- 
ing frame first law of thermodynamics for the radiation 
field, provided that one retains all the leading order terms 
with respect to P and t, including some that are of order 
P'^ , when evaluating the Lorentz transformation. 

While the equations are equivalent, the mixed frame 
formulation has two important advantages over the co- 
moving frame formulation when it comes to practical 
computation. First, we are able to write the equations in 
a manner that allows a numerical solution algorithm to 
conserve total energy to machine accuracy. We present 
such an algorithm in § |3) In contrast, it is not possible 
to write a conservative update algorithm using the co- 
moving frame equations. The reason for this is that a 
conserved total energy only exists in an inertial frame, 

^ Note that our need to retain difference between E and Eq here 
is different from the sit uatio n w hen we first apphed the Lorentz 
transformation to derive II35I I and II36I I. In that case we did not need 
to r etain the distinction between E and Eq , because in deriving 
H35| l and H36p there were no terms involving Eo expUcitly. Instead, 
Eq appeared only implicitly, as part of the flux F, and non-leading 
order corrections to F are not of leadi ng order in any regime. In 
contrast, Eq does appear explicitly in I I45I I. 
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and for a fluid whose velocity is not a constant in space 
and time, the comoving frame is not inertial. The lack 
of a conserved energy is a serious drawback to comoving 
frame formulations. 

A second advantage of the mixed-frame formulation is 
that it is far more suited to implementation in codes with 
dynamically modified grid structures such as adaptive 
mesh refinement methods. Since the radiation energy is a 
conserved quantity, it is obvious how to refine or coarsen 
it in a conservative manner. On the other hand, there 
is no obviously correct method for refining or coarsening 
the comoving frame energy density, because it will not 
even be defined in the same reference frames before and 
after the refinement procedure. 

4. AN OPTIMIZED ALGORITHM FOR STATIC DIFFUSION 
RADIATION HYDRODYNAMICS 

4.1. Operator Splitting 

Our analysis shows that for static diffusion, the terms 
involving diffusion and emission minus absorption of ra- 
diation always dominate over those involving radiation 
work and advection. In addition, some terms are al- 
ways smaller than order f3. This suggests an opportu- 
nity for a significant algorithmic improvement over ear- 
lier approaches while still retaining order (3 accuracy in 
the solution. In a simulation, one must update terms for 
the radiation field implicitly, because otherwise stability 
requirements limit the update time step to values com- 
parable to th e light-crossing time of a cell. Standard ap- 
proaches (e.g. Turner & Stonc'200r;'Whitehousc & Bat^ 
[2004 : Whitchousc ct al. 2005; Hayes ct al. 2006) there- 
fore update all terms involving radiation implicitly ex- 
cept the advection term and the radiation force term in 
the gas momentum equation. 

However, implicit updates are computationally expen- 
sive, so the simpler the terms to be updated implicitly 
can be made, the simpler the algorithm will be to code 
and the faster it will run. Since the work and advection 
terms are non-dominant, we can produce a perfectly sta- 
ble algorithm without treating them implicitly. Even if 
this treatment introduces numerically unstable modes in 
the work or advection terms, they will not grow because 
the radiation diffusion and emission/absorption terms, 
which are far larger, will smooth them away each time 
step. 

For the case of static diffusion, we therefore adopt the 
order v /c equations ^ and (|40p for mass and momentum 
conservation. For our energy equations, we adopt (j4ip 
and (|42p. but drop terms that are smaller than order (3 
for static diffusion. This gives 



-V[{pe + P)v] - Kop(47rB - cE) 



+ A 2^-1 

V«OR / 



dt 



VE 



Kop(47rS - cE) 



(47) 



A ( 2^ - 1 I V Vii; 



„ ,3 — i?2 „ 



(48) 



To solve these, we operator split the diffusion and emis- 



sion/absorption terms, which we treat implicitly, from 
the work and advection terms, which we treat explicitly. 
To do this, we write our gas/radiation state as 

i)' 

and our evolution equations as 

^ _ f , f , f . 

r^, ^c — nr I Ac— rad i ^i — radi 
Ot 

where we have broken our right hand side up into non- 
radiative terms to be handled explicitly, 

-V • (pv) 
^ _ ( -V • (pvv) - VP 



(50) 



-V • [(pe + P)v] 




(51) 



radiative terms to be handled explicitly, 

fo-rad = A(2^ - l)v • VE 



and radiative terms that must be handled implicitly, 
/ X 



(52) 



fi — rad ^ 







(53) 



-Kop(47rP - cE) 
[^■{^'^E)+noA^7rB-cE) ) 

4.2. Update Scheme 

For each update cycle, we start with the state q" at the 
old time. We first perform an implicit update to the radi- 
ation and gas energy densities using fi-rad- Any number 
of methods are possible for this. For our implementation 
of this algorithm in the Orion ada ptive mesh refinement 
AMR ) code, we use the method of iHowell fc Greenoughl 
200l . which we will not discuss in detail here. To sum- 
marize, the algorithm involves writing the equations us- 
ing second order accurate spatial discretization and a 
time discretization that limits to backwards Euler for 
large values of dE/dt (to guarantee stability) and to 
Crank-Nicolson when dE/dt is small (to achieve second 
order time accuracy). This yields a matrix equation for 
the radiation and gas energy densities at the new time, 
which may be solved on both individual grids and over a 
hierarchy of nested grids (as is necessary for AMR) us- 
ing standard multigrid techniques. The output of this 
procedure is an intermediate state q"'* which has been 
updated for fj^rad- 

Once the implicit update is done, we compute the or- 
dinary hydrodynamic update. As with the implicit up- 
date, this may be done using the hydrodynamics method 
of one's choice. For our im plementation, we use the Go- 
dunov method described bv lTruelove et al.l (|1998l ). [Kleinl 
(|1999l ). and|Fisher| ()2002[ ). This update gives us q"'t, the 
state updated for fi_rad and fe-nr- The only modification 
we make to the standard update algorithm is to include a 
radiation pressure term in the effective sound speed used 
to compute the Courant condition. Thus, we take 



Ceff 



/7P+(4/9)P(l-e 



(54) 
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and set the time step to 
At = C- 



Aa 



<;(|v| +Coff)' 



(55) 



where 7 is the ratio of specific heats for the gas, C is the 
Courant factor (usually 0.5), and the maximum is eval- 
uated over all cells. For AMR, this condition is applied 
independently on each level I, and the time step is set us- 
ing the values of At' in the standard AMR manner (e.g. 
|Klein"l999). The factor (1 — e'^oRAx^ gives us a means of 
interpolating between optically thick cells, where radia- 
tion pressure contributes to the restoring force and thus 
increases the effective signal speed, and optically thin 
cells, where radiation does not provide any pressure. 

Finally, we compute the force and advection terms in 
fc-rad- In our implementation we compute all of these at 
cell centers using second order centered differences. For 
V-B this is 



V 



2Az 



(56) 



Other derivatives are computed in an analogous manner. 
We then find the new state by 



fo-iadAt. 



(57) 



This update is manifestly only first order-accurate in 
time for the explicit radiation terms, but there is no point 
in using a more complex update because our operator 
splitting of some of the radiation terms means that we 
are performing our explicit update using a time-advanced 
radiation field, rathe r than the field at a half time step. 
(|Truelove et al.|[l998l show that one can avoid this prob- 
lem for gravitational body forces because the potential 
is linear in the density, so it is possible to derive the 
half-time step potential from the whole time step states. 
No such fortuitous coincidence occurs for the radiation 
field.) This necessarily limits us to first order accuracy 
in time for the terms we treat explicitly. However, since 
these terms are always small compared to the dominant 
radiation terms, the overall scheme should still be closer 
to second order than first order in accuracy. 

4.3. Advantages and Limitations of the Method 

Our algorithm has two significant advantages in 
comparison to other approaches, in particular those 
based on comoving frarne formulations of the equa- 
tions fe.g.lTurner fc Sto^l2001l:IWhitehouse et al.ll2005t 
iHaves et al.ll2006[ ). In any of these approaches, since the 
radiation work terms are included in the implicit up- 
date, one must solve an implicit quartic equation aris- 
ing from the combination of the terms kop{At:B — cE) 
and V: Vv. This may be done either at the same 
time one is iterating to update the flux divergence term 
V • F fWhitchousc et al. 2005), or in a separate iter- 
ation to be done once the itera tive solve for the flux 
divergence updat e is complete (jTurner fc Stond l2001l : 
iHayes et al.ll2006l ). In contrast, since our iterative up- 
date involves only K opjAnB — cE) and V • F, using the 
iHowell fc Greenoug h (2003) algorithm we may linearize 
the equations and never need to solve a quartic, lead- 
ing to a simpler update algorith m and a faster itera- 
tion step. Moreover, by using the iHowell fc Greenoughl 



(|2003f) time-centering, we obtain second order accuracy 
in time whenever E is changing slowly, as oppose d to th e 
backwards Euler d ifferencing of [Turner &: Stont (|2001h . 
Whit ehouse et all (|2 00 5X and IHaves et all (|2006f ). which 
is always first order-accurate in time. Thus, our algo- 
rithm provides a faster and simpler approach than the 
standard one. 

A second advantage of our update scheme is that it re- 
tains the total energy-conserving character of the under- 
lying equations. In each of the update steps involving ra- 
diation, for fc-rad and fi, the non-advective update terms 
in the radiation and gas energy equations are equal and 
opposite. Thus, it is trivial to write the update scheme so 
that it conserves total energy to machine precision. This 
property is particularly important for turbulent flows 
with large radiation energy gradients, s uch as those that 
occur in massive star formation (e.g. iKrumholz et al.l 
l2007f) . because numerical non-conservation is likely to 
be exacerbated by the presence of these features. In 
contr ast, in cornqving frame formalisms such as those of 
Turner fc Stond (|2001D . , Whitehouse & Bate (2004), and 
iHaves et all ^^2006^ the exchange terms in their gas and 
radiation energy equations are not symmetric. As a re- 
sult, their update schemes do not conserve total energy 
exactly. The underlying physical reason for this asym- 
metry is that total energy is conserved only in inertial 
frames such as the lab frame; it is not conserved in the 
non-inertial comoving frame. For this reason, there is no 
easy way to write a conservative update scheme from a 
comoving formulation. 

Our algorithm also has two significant limitations, one 
obvious and one subtle. The obvious limit is that our al- 
gorithm is only applicable for static diffusion problems. 
For dynamic diffusion problems, e.g. stellar interiors or 
radiation-dominated shocks, our scheme is unstable un- 
less an appropriately small timestep is used. Whether 
this instability is due to the explicit advection term, the 
explicit w ork term, or both i s not clear. Since codes such 
as ZEUS ([Haves et al.ll2006l ) treat the advection explic- 
itly without instability, however, it seems likely that the 
work term is the culprit. Regardless of the cause, even if 
we were to use a time step small enough to guarantee sta- 
bility, since the work and advection terms can be compa- 
rable to or larger than the diffusion and heating/cooling 
terms for dynamic diffusion, an algorithm that treats all 
the terms implicitly or all explicitly, rather than our mix, 
is likely to be more accurate. 

The subtle limitation is in our treatment of the hydro- 
dynamics. We perform the hydrodynamic update using a 
Riemann solver unmodified for the presence of radiation 
force, work, and heating and cooling terms. These terms 
should change the characteristic velocities of the wave 
families in ways that depend on the radiation hydrody- 
namic regime of the system. For example, in optically 
thick systems we should have a radiation-acoustic mode 
rather than a simple sound wave, and in optically thin 
systems where the radiation time scale is short compared 
to the mechanical time scale, a gas may act as if it were 
isothermal even if it has 7 1. In some cases, failure to 
modify the Riemann solver appropriately for these effects 
may produce substantial errors, including a reduction in 
the order of accuracy of the method fr om second to first 
(iPemberi Il993l : iLowrie fc More l 200 li : iMiniati fc Coiellal 
I2OO6I ). The severity of these effects for a given problem 
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depends the degree of stiffness of the radiation source 
terms. It should also be noted that the other radia- 
tion diffusion methods most commonly used for three- 
dimensional problems also suffer from this defect, so this 
is not a comparative disadvantage of our method relative 
to others. 

5. TESTS OF THE STATIC DIFFUSION ALGORITHM 

Here we describe five tests of our static diffusion algo- 
rithm, done using our implementation of the algorithm 
in the Orion AMR co de, various aspects of whic h are 
described in detail by iPuckett fc Saltznianl (Il992l mul- 
tifluid hydrodynamics), 'Truelove et al.' ("lOQS', hydrody- 
namics and gravity), Klein (1999, hydrodynam ics and 
gravit y), iFished (j2002L gr avitvl. IHowell fc Greenoughl 
(|2003l radiatio n transport). iKrumholz et al.l ( 2004 sink 
particles), and iCrockett et al.l ()2005l . magnetohvdrodv- 
namics). For all of these tests we use a single fluid with 
no magnetic fields and no self-gravity. 

5.1. Non- Equilibrium Marshak Wave 

As an initial check of the gas-radiation energy exchange 
in our code in a case when radiation pressure is not 
significant and the gas is at rest, we simulate the non- 
equilibrium Marshak wave problem. In this problem, a 
zero-temperature, motionless, gaseous medium occupy- 
ing all space at z > is subject to a constant radiation 
flux i<]ncZ incident on its surface at z = 0. The gas is 
held stationary, appropriate for early times before hy- 
drodynamic motions become significant. The medium is 
gray, with opacity kor = kqp = k, and the constant- 
volume specific heat capacity of the gas is taken to have 
the same dependence as that of the radiation, i.e. 
c„ = [d{e - v^/2)/dTg]^ = aT^, where Tg is the gas 
temperature. The gas is not assumed to be in thermal 
equilibrium with the radiation field, so the gas and radi- 
ation te mpera t ures m ay be different. 

ISu fc OlsonI (|1996f ) give a semi-analytic solution to the 
time-dependent behavior of the radiation energy density 
E{z,t) and gas temperature Tg{z,t) for this problem. 
They introduce the dimensionless position and time vari- 
ables X = -\/3kz and r = {AaRcn/ a)t, and the "retarda- 
tion" parameter e = Aan/a, and show that the dimen- 
sionless radiation energy density 



and 



w(x,r)=(^) 



E{z,t) 



(58) 



= 1 
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The dimensionless gas energy density is 
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(65) 



Numerical evaluation of the integrals ([55]) and (|64p for 
u and V is not trivial because the integrands perform 
an infinite number of oscillations about zero as 77 ^ 0. 
Correct computation of the result when r is small and x 
is large requires careful numerical analysis to ensure that 
the positive and negative contributions cancel properly 
(J. Bolstad, 2007, in preparation). 

We compare the properly computed semi-analytic re- 
sults for u and w to a calculation performed with Orion 
using K = 1 cm~^ and a = 32ar/c (so e = 0.5). The 
computational domain goes from to 15 cm (and thus 
to an optical depth kz = 15), and is resolved by 100 
equally-sized cells. For this test, since we are comparing 
to a pure diffusion result, we set the flux limiter A = 1/3 
everywhere. 

Figures [1] and [2] compare the semi-analytic dimension- 
less radiation and gas energy densities with the values 
computed by Orion. At r = 0.001 the agreement is fairly 
poor due to low numerical resolution, since the wave only 
reaches an optical depth of kz ~ 0.2 and kz = 0.15 is the 
size of an individual computational cell. However, at 
later times when the wave is resolved by a reasonable 
number of cells, the agreement between the code result 
and the semi-analytic solution is excellent. 

5.2. Radiating Blast Wave 

We next compare to a test problem in which the gas is 
not at rest: a Sedov-type blast wave with radiation dif- 
fusion. iReinicke fc Mever-ter-Vehiil (|1991h gave the first 
similarity solution to the problem of a point explosion 
wit h heat conduction, and f ollowing Shestakov (199!|) 
and IShestakov fc Greenoughl (j200lD . we can adapt this 
solution to the case of a point explosion with radiation 
diffusion. This tests our code's ability to follow coupled 
radiation-hydrodynamics in cases where radiation pres- 
sure is small. 

We first summarize the semi-analytic solution. Con- 
sider an 7i = 3 dimensional space filled with an adiabatic 
gas with equation of state P — {"/ ~ l)/3e = ^pT, where 
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Fig. 1. — Dimensionless radiation energy density u versus optical 
depth Kz at a series of times r. We show the semi-analytic solution 
(solid lines) and the solution computed with Orion {dashed lines). 
The value of the dimensionless time r is indicated by each curve. 
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Fig. 2. — Same as Figure[T] but for the dimensionless gas energy 
density v. 

r is the gas constant. The Planck mean opacity kqp of 
the gas is very high, so the gas and radiation temper- 
atures are always equal. The Rosseland mean opacity 
has a powerlaw form kqr = i^oR,oP"^''T~"', and we assume 
that it is always high enough to place us in the diffusion 
regime, so A = 1/3. Note that the choice of —n = —3 
as the exponent of the opac ity powerlaw is a necessary 
condit ion for applying the iReinicke fc Mever-ter-Vehnl 
(|1991| ) conduction solution to our radiation diffusion 
problem. Moreover, the similarity solution does not in- 
clude radiation energy density or pressure, so we consider 
only temperatures for which the gas energy density and 
pressure greatly exceed the radiation energy density and 
pressure, i.e. pe ^ (irT*. 

Under the assumptions described above, we may re- 
write the gas and radiation energy equations (|4ip and 
([1^ as a single conduction-type equation for the tem- 
perature, 

-T = V(xop'^r^VT), (66) 



dt 



where = de/dT = r/(7 — 1) is the constant- 
volume specific heat of the gas, xo = 4cafl'/(3Kofl,o), 
a — —m, and b = n + 3. This equation has the 



same form as the conduct i on eq uation considered by 
IReinicke fc Mever-ter-Veh^ (|1991| ). 

Consider now a point explosion at the origin of a 
spherically symmetric region with an initial powerlaw 
density distribution p{r,t = 0) = gor~'''' . Initially 
the gas temperature T and pressure P are negligible. 
The explosion occurs at the origin at time zero, so the 
initial gas energy density i s (pe) (r,t = 0) = Eo6{r). 
IReinicke fc Mever-ter-Vehnl (|1991t ) show that if the ini- 
tial density profile has a powerlaw index 



_ (26 - l)n + 2 
26-2a+l ' 



(67) 



then one may obtain a similarity solution via the change 
of variables 



G(e) = 



r 

P{r,t) 



t 



ar 



e(e)=T(r,t)r( 



t) 



(68) 
(69) 
(70) 
(71) 



Here, ^, G(^), U{(,), and 6(^) are the dimensionless dis- 
tance, density, velocity, and temperature, 



26 - 2a -f 1 
26- {n + 2)a + n 



(72) 



and C is a constant with units of [length] [time] whose 
value is determined by a procedure we discuss below. 

With this similarity transformation, the equations of 
motion and heat conduction reduce to 



[/' - (1 - C/)(lnG)' + (n - fcp)[/: 
(1 - U)U' + U{a-'^ - J7): 







(73) 



(74) 



and 



(75) 



2[U' + nU- ^l{a-^ - 1)] = ^^(1 - U)[\u{^'^ Q)]' 

+ poO'^G'^-^c'-^''-^^^'" ■ ((ine)" + M^^e)]' 

■{n-2 + a[ln(r'^-G)]' + (6 + l)[ln(e2e)]'}) ,(76) 
where ()' = d{)/d\n^, ^ = 2/(7 - 1), and 
2xo(aCi/")2^- 



-sgn(<). 



(77) 



This constitutes a fourth-order system of non-linear or- 
dinary differential equations. All physical solutions to 
these equations pass through two discontinuities, a heat 
front and a shock front, with the heat front at larger 
radius. However, the jump conditions for these discon- 
tinuities are easy to determine, and one can integrate 
between them. For a given /3o, the solution depends only 
on the dimensionless parameter 



n = 



2X0 



Eo 
90 



6-1/2 



(78) 



which measures the strength of the explosion. Large val- 
ues of constitute "strong" explosions, and the ratio of 
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heat front radius to shock front radius is a monotonicaUy 
increasing function of fl. It is important at this point to 
add a cautionary note: in deriving the similarity solution, 
we assumed that radiation energy density is negligible in 
comparison to gas energy density. This cannot strictly 
be true at early times, since at i = the temperature 
diverges at the origin, and the radiation energy density 
varies as T to a higher power than the gas energy den- 
sity. However, the true behavior should approach the 
similarity solution at later times. 

While we have reduced the gas dynamical equations to 
a system of ordinary differential equations that is trivial 
to integrate, solving the full problem is complex because 
the equations still depend on the unknown parameter 
Po, which in turn depends on C. To solve the prob- 
lem, we must determine Pn from the give n initial con- 
ditions. iReinicke fc Mever-ter-VehnI (|199lD describe the 
iteration procedure required to do this in detail, and we 
only summarize it here. To find a solution, one first 
chooses a value > 1 for the dimensionless radius of 
the heat front, applies the boundary conditions at the 
front, and guesses a corresponding value of /3o- For each 
£,h there exists a unique /3o for which it is possible to 
integrate the equations back from £, — £,h the loca- 
tion of the shock front at ^ = ^s, apply the shock jump 
conditions, and continue integrating back to the origin 
at ^ = without having the solution become double- 
valued and thus unphysical. One iterates to identify 
the allowed value of /3o for the chosen ^h, and this gives 
the unique density, velocity, and temperature profiles al- 
lowed for that S,h- However, the solution one finds in 
this way may not correspon d to the desired value of fl. 
IReinicke fc Mever-ter-Vehnl show that 



27r 



6-1/2 



(79) 



Thus, each choice of corresponds to a particular value 
of fl, and one must iterate a second time to find the value 
of that gives the value of Q determined from the input 
physical parameters of the problem. Alternately, instead 
of specifying a desired value of 17, one may specify a ratio 
R = £.h/£,s, which also determines a unique value for ■ 
For our comparison between the semi-analytic solu- 
tion and Orion, we adopt the parameters 7 = 7/5, 

Cv = 1/(7 -!),« = -2, 6 = 6, go = Xo = 1, and 
Eq = 135, which yields a strength fl = 1.042 x 10^^ and 
a ratio R = 2.16. In the simulation, we turn off terms 
in the code involving radiation pressure and forces, and 
we set A = 1/3 exactly. We use one-dimensional spheri- 
cal polar coordinates rather than Cartesian coordinates; 
the solution procedures for this are identical to the ones 
outlined in § 21 with the exception that the gradient 
and divergence operators have their spherical rather than 
Cartesian forms, and the cell-centered finite differences 
are modified appropriately. Our computational domain 
goes covers < r < 1.05, resolved by 256, 512, or 1024 
cells, and has reflecting inner and outer boundary condi- 
tions. To initialize the problem we set initial density to 
the powerlaw profile p = r~^p (with kp set from equa- 
tionl67p. the initial velocity to zero, and the initial energy 
density to a small value, except in the cell adjacent to 
the origin, where its value is pe = 135/(7 ~ !)• 
Figures [3l IH and [5] compare the semi-analytic density. 
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Fig. 3. — Density p versus radius r for the radiating blast wave 
test. We show the semi-analytic solution {solid line), and the Orion 
results at resolutions of 256, 512, and 1024 cells {dashed lines). 
The 256-cell run is the dashed line furthest from the semi-analytic 
solution, and the 1024-cell run is the dashed line closest to it. 
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Fig. 4. — Same as Figure [S] but for the velocity v. 

velocity, and temperature profiles to the values we obtain 
from Orion after running to a time t = 0.06. As the 
plots show, the Orion results agree very well with the 
semi-analytic solution, and the agreement improves with 
increasing resolution. In the lowest resolution run, there 
is a small oscillation in the density and velocity about a 
third of the way to the shock, which is likely due to the 
initial blast energy being deposited in a finite-volume 
region rather than as a true S function. However, this 
vanishes at higher resolutions. Overall, the largest errors 
are in the temperature in the shocked gas. 

As a metric of convergence, we plot the error of our 
simulation relative to the analytic solution as a function 
of resolution in Figure [51 We do this for the quantities 
r/i and Tj, the positions of the shock and heat fronts, and 
their ratio R. For this purpose, we define the location of 
the heat and shock fronts for the simulations as the po- 
sitions of the cell edges where dT/dr and dp/dr are most 
negative. As the plot shows, at the highest resolution the 
errors in all three quantities are ^ 3%, and the calcula- 
tion appears to be converging. The order of convergence 
is roughly 0.6 in all three quantities. It is worth not- 
ing that computing the locations of the heat and shock 
fronts is a particularly strong code test, because obtain- 
ing the correct propagation velocities for the two fronts 
requires that the code conserve total energy very well. 



14 



Krumholz, McKee, & Klein 




Fig. 5. — Same as Figure [3l but for the temperature T. 




Fig. 6. — Fractional error versus resolution N in the ra- 
diating blast wave test. The fractional error is defined as 
simulation value — analytic value | /analytic value. We show error 
in the heat front radius (plus signs), shock front radius (as- 
terisks), and their ratio R = r^/rs (diamonds) . 



Non-cons ervative codes have significant difRcultics with 
this test (|Timmes et aLll2006D . 

5.3. Radiation Pressure Tube 

Our third test is to simulate a tube filled with radia- 
tion and gas. The gas within the tube is optically thick, 
so the diffusion approximation applies. The two ends of 
the tube are held at fixed radiation and gas temperature, 
and radiation diffuses through the gas from one end of 
the tube to the other. The radiation flowing through 
the tube exerts a force on the gas, and the gas density 
profile is such that, with radiation pressure, the gas is in 
pressure balance and should be stationary. For computa- 
tional simplicity, we set the Rosseland- and Planck-mean 
opacities per unit mass of the gas to a constant value k. 
A simulation of this system tests our code's ability to 
compute accurately the radiation pressure force in the 
very optically thick limit. 

We first derive a semi-analytic solution for the con- 
figuration of the tube satisfying our desired conditions. 
Since the gas is very optically thick and we are starting 
the system in equilibrium, we set Tiad = Tgas = T. The 
fluid is initially at rest. The condition of pressure balance 
amounts to setting d{pw)/ dt 4- V • (pvv) = in equation 
pO|) , so that the radiation pressure force balances the gas 
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Fig. 7. — Density, temperature, and pressure versus position in 
the radiation tube problem. The bottom panel shows total pressure 
(solid line), gas pressure (dashed line), and radiation pressure (dot- 
dashed line). 

pressure gradient. Thus, we have 

dP 



dx 



+ A 



dE 
dx 



= 



dx /i da; 



(80) 
(81) 



In the second step we have set E = aaT'^ and P = 
pkBT/ fi, where /i is the mean particle mass, and we have 
set A = 1/3 as is appropriate for the optically thick limit. 
The radiation energy equation (|48p for our configuration 
is simply 



dx^ 



V [dx 



d 

dx 

1 rdp 

p \dx 



cA dE\ _^ 
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dx 



(82) 



(83) 



Equations ([5T|l and 



are a pair of coupled non- 
linear ordinary differential equations for T and p. The 
combined degree of the system is three, so we need 
three initial conditions to solve them. Thus, let the 
tube run from x = xq to x = xi, with temperature, 
density, and density gradient Tq, po, and {dp/dx)o at 
xq. For a given choice of initial conditions, it is triv- 
ial to solve ([5T|) and ([55)1 numerically to find the den- 
sity and temperature profile. We wish to investigate 
both the radiation pressure and gas pressure dominated 
regimes, so we choose parameters to ensure that our 
problem covers both. The choice a;o = 0, a:i = 128 
cm, pq — 1 g cm^^, {dp/dx)o — b x 10^^ g cm^^, and 
To = 2.75 X 10^ K satisfies this requirement if we adopt 
/i = 2.33 TOP = 3.9 X 10~24 g and K = 100 cm^ g-i. Fig- 
ure [7] shows the density, temperature, and pressure as a 
function of position for these parameters. 

We solve the equations to obtain the density and tem- 
perature as a function of position, and then set these val- 
ues as initial conditions in a simulation. The simulation 
has 128 cells along the length of the tube on the coarsest 
level. We impose Dirichlet boundary conditions on the 
radiation field, with the radiation temperature at each 
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Fig. 8. — Relative error in density (solid line), gas temperature 
(dashed line), and radiation temperature (dot-dashed line) in the 
radiation tube test. 

end of the tube set equal to its value as determined from 
the analytic solution. We use symmetry boundary con- 
ditions on the hydrodynamics, so that gas can neither 
enter nor leave the computational domain. To ensure 
that our algorithm does not encounter problems at the 
boundaries between AMR levels, we refine the central 
1/4 of the problem domain to double the resolution of 
the base grid. We evolve the system for 10 sound cross- 
ing times and measure the amount by which the density 
and temperature change relative to the exact solution. 
We plot the relative error, defined as (numerical solution 
— analytic solution) / (analytic solution), in the density, 
gas temperature, and radiation temperature in Figure [51 
As the plot shows, our numerical solution agrees with 
the analytic result to better than 0.5% throughout the 
computational domain. The density error is smallest in 
the higher resolution central region, as expected. There 
is a very small increase in error at level boundaries, but 
it is still at the less than 0.5% level. 

5.4. Radiation- Inhibited Bondi Accretion 

The previous test focuses on radiation pressure forces 
in the optically thick limit. To test the optically thin 
limit, we simulate accretion onto a radiating point par- 
ticle. We consider a point mass M radiating with a con- 
stant luminosity L accreting from a background medium. 
The medium consists of gas which has zero velocity and 
density poo far from the particle. We take the gas to 
be isothermal with constant temperature T, and enforce 
that it is not heated or cooled radiatively by setting its 
Planck opacity kqp = 0. We set the Rosseland opacity of 
the gas to a constant non-zero value kor, and choose p^o 
such that the computational domain is optically thin. 
In this case, the radiation free-streams away from the 
point mass, and the radiation energy density and radia- 
tive force per unit mass on the gas are 



E- 



L 



A'Kr'^c 
kqrL 



(84) 
(85) 



where r is the radial vector from the particle and r is 
its magnitude. The gravitational force per unit mass is 
fg = —{GM/r'^){r /r), so the net force per unit mass is 

f = f. + fg = -(l-/Edd)^0, (86) 
where ^ 

is the fraction of the Eddington luminosity with which 
the point mass is radiating. 

Since the addition of radiation does not alter the 1/r^ 
dependence of the specific force, the solution is simply 
the standard Bondi (1952) solution, but for an effective 
mass of (1 — /Edd)Af- The accretion rate is the Bondi 
rate 



Mb = 4<r^c,p„ 



(88) 



where 



(1-/1 



GM 



EddJ 



is the Bondi radius for the effective mass, Cg is the gas 
sound speed at infinity, and ^ is a numerical factor of 
order unity that depends on the gas equation of state. 
For an isothermal gas, ^ = e'^/^/4, and the radial profiles 
of the non-dimensional density a = p/poo and velocity 
u = v/cs are given by the so lutions to the non-linear 
algebraic equations (|Shulll99^ 



x^au : 
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(90) 
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where x = r/rs is the dimensionless radius. 

To set up this test, we make us e of the Lagrangian sink 
particle algorithm of fKrumholz et al.l (I2004D. coupled 
with the "star particle" algorithm of iKrumholz et al.l 
()2007f l which allows the sink particle to act as a source 
of radiation. We refer readers to those papers for details 
on the sink and star particle algorithms. We simulate a 
computational domain 5 x 10^^ cm on a side, resolved 
by 256'^ cells, with a particle of mass M = 10 Mq and 
luminosity L — 1.6 x 10^ Lq at its center. We adopt 
fluid properties poo — 10^^* g cm~'^, kqr = 0.4 cm^ 
g~^, and Cs = 1.3 x lO'' cm s~^, corresponding to a gas 
of pure, ionized hydrogen with a temperature of 10^ K. 
With these values, /Edd = 0.5, = 4.0 x 10^^ cm, and 
Mb = 2.9 X 10"'^'^ g s~^. We use inflow boundary con- 
ditions on the gas and Dirichlet boundary conditions on 
the radiation field, with the radiation energy density on 
the boundary set to the value given by equation (j84p . 

Figure [5] compares the steady-state density a and ve- 
locity u computed by Orion to the analytic solution. The 
agreement is excellent, with differences between the an- 
alytic and numerical solutions of ~ 1% everywhere ex- 
cept very near the accretion radius at a; = 0.25. The 
maximum error is ~ 10% at the surface of the accre- 
tion region; this is comparable to the error in density 
for no n-radiative Bond i accre tion with similar resolu- 
tion in IKrumholz eFall ()2004D . In comparison, the so- 
lution is nowhere near the solution that would be ob- 
tained without radiation. After running for 5 Bondi 
times {— rs/cs), the average accretion rate is 2.4 x 10^^ g 
s~^. While this differs from the analytic solution by 19%, 
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Fig. 9. — Dimensionless density a (upper panel) and veloc- 
ity u {lower panel) versus dimensionless position x for radiation- 
inhibited Bondi accretion. We show the analytic solution {solid 
line), the solution as computed with Orion {dashed line), and 
the analytic solution for Bondi accretion without radiation {dotted 
line). For the Orion result, the values shown are the radial aver- 
ages computed in 128 logarithmically-spaced bins running from the 
accretion radius x = 0.25 to the outer edge of the computational 
domain a; = 5. 

the error is also not tremendously different from that ob- 
tained bv lKrumholz et al] (|2004l ) when the Bondi radius 
was resolved by 4 accretion radii, and is nowhere near 
the value of 1.2 x 10^* g s~^ which would occur without 
radiation. 

We should at this point mention one limitation of our 
algorithm, as applied on an adaptive grid, that this test 
reveals. The l/r^ gradient in the radiation energy den- 
sity is very steep, and we compute the radiation force by 
computing gradients in E. We found that, in an AMR 
calculation, differencing this steep gradient across level 
boundaries introduced significant artifacts in the radia- 
tion pressure force. With such a steep gradient, we were 
only able to compute the radiation pressure force accu- 
rately on fixed grids, not adaptive grids. This is not a sig- 
nificant limitation for most applications though, since for 
any appreciable optical depth the gradient will be much 
shallower than 1/r^. As the radiation pressure tube test 
in 15.31 demonstrates, in an optically thick problem the 
errors that arise from differencing across level boundaries 
are less than 1%. 

5.5. Advecting Radiation Pulse 

The previous two tests check our ability to compute 
the radiation pressure force accurately in the optically 
thick and optically thin limits. However, they do not 
strongly test radiation advection by gas. To check this, 
we simulate a diffusing, advecting radiation pulse. The 
initial condition is a uniform background of gas and ra- 
diation far from the pulse. Centered on a:: = there is 
an increase in the radiation energy density and a corre- 
sponding decrease in the gas density, so that the initial 
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Fig. 10. — Density, temperature, and pressure versus position 
in the advected radiation pulse problem. The bottom panel shows 
total pressure (solid line), gas pressure (dashed line), and radiation 
pressure (dot-dashed line). 

condition is everywhere in pressure balance. As radiation 
diffuses out of the pulse, pressure support is lost and the 
gas moves into the lower density region. We cannot solve 
this problem analytically, but we can still perform a very 
useful test of the methodology by comparing a case in 
which the gas is initially at rest with respect to the com- 
putational grid with a case in which the gas is moving at 
a constant velocity with respect to the grid. The results 
should be identical when shifted to lie on top of one an- 
other, but the work and advection terms will be different 
in the stationary case than in the advected case. Check- 
ing that the results do not change when we advect the 
problem enables us to determine if our code is correctly 
handling the advection of radiation by the gas. 

For our simulations, we use equal initial gas and radia- 
tion temperatures, with temperature and density profiles 

r = To + (Ti-To)exp(^-^^ (92) 

with To = 10^ K, Ti = 2 X 10^ K, po = 1.2 g cm^^^ w = 
24 cm, fi = 2.33 TOP = 3.9 x 10~24 and k = 100 cm^ 
g^^. The density, temperature, and pressure profiles are 
shown in Figure [TUl In the bottom panel, the solid line 
is the total pressure, the dashed line is the gas pressure, 
and the dot-dashed line is the radiation pressure. As 
the figure indicates, the system is initially in pressure 
balance. 

We compare two runs, one where the velocity is zero 
everywhere and another with a uniform initial velocity 
V = 10® cm s~^ in the x direction. In both runs the sim- 
ulation domain extends from —512 to 512 cm, resolved 
by 512 cells with no adaptivity. We use periodic bound- 
ary conditions on the gas and the radiation, and run for 
4.8 X lO^'^ s, long enough for the pulse to have been ad- 
vected over its own initial width twice. 

To check our results, we shift the advected run by 
48 cm in the —x direction, so that it should lie on top 
of the unadvected run. Figure [TT] shows the configura- 
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Fig. 11. — Density, velocity, and temperature in the adverted 
radiation pulse problem after 4.8 X 10~ s of evolution. In each 
panel the solid line is the unadvected run, and the dashed line is 
the adverted run shifted 48 cm in the —x direction. In the velocity 
plot, the velocity we show for the adverted run is relative to the 
overall systematic velocity of 10^ cm in the initial condition. 



0.02 




-0.02 r , , , I , I ... I ... I ... I ... ' 

-600 -400 -200 200 400 600 

X (cm) 

Fig. 12. — Relative error in density (solid line) and gas/radiation 
temperature (dashed line) in the radiation pulse test. 

tion of the advected and unadvected runs at this point. 
We then plot the relative difference between the ad- 
vected and unadvected runs, defined as (unadvected — 
advected) /unadvected, in Figure [121 We do not differ- 
entiate between the gas and radiation temperatures, be- 
cause they are identical at the 10"'^ level. We do not 
plot the error in velocity because the velocities in the 
unadvected run are close to zero over most of the com- 
putational domain. As the plot shows, the difference 
between the advected and unadvected runs is less than 
2% everywhere in the simulation. 

6. SUMMARY 

We derive the correct equations for mixed frame flux- 
limited diffusion radiation hydrodynamics. The error in 



our equations if of order j (? in the static diffusion limit, 
and of order v j c in the dynamic diffusion and streaming 
limits. We give the equations in a form that is well- 
suited to implementation in numerical simulations, be- 
cause they make it trivial to maintain exact conservation 
of total energy. Our analysis reveals that lower order for- 
mulations of the equations, which neglect differences be- 
tween the laboratory and comoving frames, are incorrect 
at order unity for systems in the dynamic diffusion limit. 
It remains to be seen how serious this defect is in prac- 
tice, but analytic arguments suggest that at a minimum 
one ought to be very careful in applying zeroth order 
codes to problems where there are interesting or impor- 
tant structures on scales for which /?t ~ 1. We give the 
equations that are correct to leading order for dynamic 
diffusion, which do not suffer from this problem. 

Our analysis also reveals that, for static diffusion prob- 
lems, one can obtain a significant algorithmic simplifica- 
tion and speedup compared to algorithms based on co- 
moving frame formulations of the equations by treating 
non-dominant radiation terms explicitly rather than im- 
plicitly. This advance is possible even though the under- 
lying equations of our method conserve total energy to 
machine precision while comoving frame formulations of 
the equations do not. This property is particularly im- 
portant for flows that are turbulent or otherwise involve 
large gradients in gas or radiation properties, since these 
are the problems most likely to suffer from numerical 
non-conservation. We demonstrate an implementation of 
this method in the Orion adaptive mesh refinement code, 
and show that it provides excellent agreement with ana- 
lytic and semi-analytic solutions in a series of test prob- 
lems covering a wide range of radiation-hydrodynamic 
regimes. 
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APPENDIX 



A. SCALINGS IN THE DYNAMIC DIFFUSION LIMIT 

Here we show that the emission minus absorption term AttB/c — _E is of order P'^E in the dynamic diffusion hmit. 
IMihalas k Weibel-Mihalal (fl999l) argue that in this hmit AttB/c — £' is of ordei (P/t)E. However, this conclusio n is 
based on their analysis of the second order equilibrium diffusion approximation (IMihalas fc Weibel -Mihalasl ll999l pg. 
461-466), in which they retain terms of order /3/r while dropping those of order . While this is correct for static 
diffusion, in the dynamic diffusion limit ^ P/t, so the approach in IMihalas fc Weibel-MihalasI is not consistent, 
and is insensitive to terms of order 

We will not give a general proof that 4ttB/c ~ E ^ P'^E ior dynamic diffusion, but we can establish it by a simple 
thought experiment. Consider a system that is infinitely far into the dynamic diffusion limit, in the sense that t = oo: 
an infinite uniform medium that is at rest and in perfect thermal equilibrium between the radiation field and the gas. 
In the rest frame of the medium, these assumption require E'o = 4:TtB/c, Fq = 0, and Vq — {Eq/3)T. Now consider 
an observer moving at velocity v relative to the medium. In the observer's frame, A-kB/c is the same because the gas 
temperature Tq is a world-scalar, and the Lorentz transform to all orders for the energy gives 



E = Y \ Eo + 2 



= 7 



.vFq 

2 



(vv) 



1 + - ^ 



(AttB 



AttB^^ 



(Al) 
(A2) 
(A3) 



Thus, for this case it is clear that AttB/c — E ^ P^E to leading order. 

Note that using the correct scaling is necessary to obtain sensible behavior from the equations in the dynamic 
diffusion limit. If one assumes that AttB/c — E ^ {P/t)E, then in the gas and radiation energy equations (j4T]l and 
([1^ in the dynamic diffusion limit, the term Kop(f^/c)[(3 — R2)/2]E is of higher order than any other term except 
perhaps the time derivative. Since this term is non-zero for any system with non-zero velocity, opacity, and radiation 
energy density, this means that there would be no way for the time derivative term to ever vanish. Thus, a system in 
the dynamic diffusion limit could never be in equilibrium unless its velocity or radiation energy were zero everywhere. 
Clearly this cannot be correct, since it predicts that our static, infinite, uniform medium cannot be in equilibrium 
when seen by an observer moving by at velocity v, even though it is manifestly in equilibrium in its own rest frame. 
On the other hand, if we take AttB/c — E = (4/3)(u^/c^)i?, as computed from the Lorentz transform, it is trivial to 
verify that equations (|¥T|) and (|^^ correctly give d{pe)/dt = dE/dt — 0, and (G°, G) = (0, 0) as well. The observer 
sees a flux that does work on the gas, but this is precisely canceled by a mismatch between emission and absorption 
of radiation by the gas, leading to zero net energy transfer. 



